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ABSTRACT 

Aims. We investigated the response of the solar atmosphere to non-thermal electron beam heating using the radiative transfer and 
hydrodynamics modelling code RADYN. The temporal evolution of the parameters that describe the non-thermal electron energy 
distribution were derived from hard X-ray observations of a particular flare, and we compared the modelled and observed parameters. 
Methods. The evolution of the non-thermal electron beam parameters during the X1.5 solar flare on 2011 March 9 were obtained 
from analysis of RHESSI X-ray spectra. The RADYN flare model was allowed to evolve for 110 seconds, after which the electron 
beam heating was ended, and was then allowed to continue evolving for a further 300s. The modelled flare parameters were compared 
to the observed parameters determined from extreme-ultraviolet spectroscopy. 

Results. The model produced a hotter and denser flare loop than that observed and also cooled more rapidly, suggesting that additional 
energy input in the decay phase of the flare is required. In the explosive evaporation phase a region of high-density cool material 
propagated upward through the corona. This material underwent a rapid increase in temperature as it was unable to radiate away all 
of the energy deposited across it by the non-thermal electron beam and via thermal conduction. A narrow and high-density (iig < lO'^ 
cm“^) region at the base of the flare transition region was the source of optical line emission in the model atmosphere. The collision¬ 
stopping depth of electrons was calculated throughout the evolution of the flare, and it was found that the compression of the lower 
atmosphere may permit electrons to penetrate farther into a flaring atmosphere compared to a quiet Sun atmosphere. 
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1. Introduction 


Non-thermal electrons are thought to be accelerated in the 
corona during solar flares and lose the majority of their energy in 
the lower atmosphere through electron-electron Coulomb colli¬ 
sions. A small amount of their energy is lost through electron-ion 
bremsstrahlung radiation, producing hard X-ray (HXR) photons. 
If the accelerated electrons are thermalised in the unresolved 
HXR source region within the observational interval, then the 
emissio n is referred to as thick-target. Collisional thick-target 
models (iBrow ^ ll97lUHuds^ll972i lR row ^ Il973l) were origi¬ 
nally proposed to explain HXR bursts and have since become 
the standard model to explain energy transfer into and heating of 
the low solar atmosphere during the flare impulsive phase. 

The hydrodynamics and radiation of a loop s ubject to thick- 
target heating was studied b\ iFisher et alJ (Il985clf9la) . who found 
two distinct types of hydrodynamic response to beam heat¬ 
ing based on the non-thermal electron energy flux, referred 
to as gentle and explosive evaporation. Gentle evaporation oc¬ 
curs when the heating rate in the chromosphere is insufficient 
to increase the temperature past the peak in the radiative loss 
function, while explosive evaporation occurs when the heating 
timescale is short e r than the hydrodynamic expansion timescale. 
iHawlev & Fisher! (Il994l) included the non-local therm odynamic 
equili brium (NLTE) radiative transfer code MULTI (ICarlssonl 
Il986l) in their solar flare models to calculate optically thick tran¬ 
sitions from hydrogen, singly ionised calcium, and magnesium, 
but limitations were imposed on the hydrodynamics of the atmo¬ 


sphere. The RA DYN code (ICarlsson & Steinl[T9^[I^[l9^ 
was adopted by lAbbett & HawlevI (Il999t) to study the detailed 
radiative hydrodynamics of a loop subjected to el ectron beam 
heatin g. Improvements to the code were made by lAlIred et~an 
(I 2 OO 5 I) to include heating from an injected electron energy spec¬ 
trum characterised by a double-power-law distribution and time- 
dependent beam heating parameters. Additional transitions were 
also included in the treatment of soft X-ray (SXR), extreme- 
ultraviolet (EUV), and ultra-violet (UV) heat ing, referred to as 
radiative back-warming (iMachado et al.lT98^ . 


There are recent examp les in the lit e rature of applications 
of the code to flare studies. ICheng et al.l (l2010h employed RA¬ 
DYN to study continuum emissions in white-light solar flares 
as a function of heliographic angle and energy deposition rate 
into a sunspot atmosphere. It has also been used to compare the 
observed and pred icted intensity in emission fin es and passbands 
during solar flares. iRubio da Costa et alj(l2012l) examined the in¬ 
tensity of the Her and Lya emission fines during several flares 
and compared tfie results witfi the synthetic line intensities from 
a number of RADYN models. They found that the predicted Ha 
intensity matched with observations to within an order of mag¬ 
nitude and agr eed better wi t h the o bservations than the predicted 
Lya intensity. ICheng et alJ (l2012h employed RADYN to create 
synthetic UV fight curves by convolving the modelled UV con¬ 
tinuum with the TRACE 1600A passband. The modelled emis¬ 
sion decayed faster than the observations, though the heating 
lasted for only 2.5s and the contribution of C iv emission lines to 
the filter were not considered. Radiative hydrodynamic simula- 
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tions have applications to th e study of stellar flares d Allred et al.l 
l2006tlifowalski et al.l2013h as there is a greater wealth of optical 
spectra available for investigation compared to solar flares. 

Previous studies involving flare modelling have typically 
used general parameters that characterise the non-thermal elec¬ 
tron energy distribution. However, recent multi-wavelength ob¬ 
servations of solar flares from X-ray to optical wavelengths have 
motivated detailed radiative hydrodynamical modelling of flares 
with realistic, time-dependent electron beam heating functions 
derived from HXR observations, such a s the detailed s t udy o f 
the X2.2 event on 2011 February 15 bv iMilligan et al.l (l2014l) . 
The aims of such studies are to investigate the chromospheric 
response to electron beam heating, determine which line and 
continuum transitions dominate the radiative energy losses, and 
directly compare the model results to optical, UV, and EUV ob¬ 
servations of the same event. 

In this study we model the atmospheric response to non- 
thermal electron beam heating using RADYN with the param¬ 
eters of the injected electron energy spectrum derived from 
RHESSI HXR observations of a specific flare. The results of 
the model are then directly compared to physical parameters de¬ 
rived by other instruments, including hydrodynamic variables 
and radiative output of EUV emission lines and continua ob- 
served by the Extre me-ultraviolet Variability Experiment (EVE, 
IWoods et al.ll2012h . In Sect. we describe the analysis of the 
HXR emission used to determine the input parameters to the 
modelling code, while in Sect. [3 the calculation and evolution 
of the model atmosphere is outlined. Section |4] shows parame¬ 
ters of the flare derived from observations and compares these 
with the model calculations. A discussion of the results and our 
conclusions are presented in Sect. |5] 


2. Observations 

The event analysed as part of this study was a Geostationary 
Orbiting Environmental Satellites (GOES) XI.5 class flare that 
occurred in the Pyd NOAA active region 11166 (N08W11) on 
2011 March 9 (SOL201 l-03-09T23:23). It was observed by 
instruments onboa rd the Solar Dynamics Observatory (SDO, 
iPesnell et al.ll2012h and by the Reuv en Ramaty High Energy So¬ 
lar Spectroscopic Imager (RHESSI. fLin et aHl2002h until it en¬ 
tered the night-time portion of its orbit at approximately 23:38 
UT. Increased pre-flare SXR emission from thermal plasma was 
observed starting at 23:15 UT, and increased HXR emission up 
to energies of 100 keV were detected from 23:20 - 23:22 UT 
(see Eig. [T]i. Observations from RHESSI were employed to de¬ 
termine the non-thermal electron energy distribution and HXR 
source area of the flare, with the aim of adopting these parame¬ 
ters as input to the modelling code RADYN. 


2.1. X-ray imaging 


To derive the energy deposition rate into the lower atmosphere, 
it was necessary to estimate the footpoint area of the event. 
The area of the footpoint HXR emission was estimated using 
the imaging capabilities of RHESSI. Initial images of the ther¬ 
mal and non-thermal emissio ns were made using the Clean al¬ 
gorithm dHurford et al.ll20()3) over a number of energy bands. 
However, for esti mates of the HXR f ootpoint source area, the 
Pixon algorithm (iMetcalf et al.l 199^ was used. In a study of 
HXR source sizes iDennis & PernakI (1200^ found that the Pixon 
algorithm, the maximum-entropy method, and Clean compo¬ 
nents gave similar estimates of the HXR area. A similar study 



Start Time (09-Mar-11 23:00:12) 

Fig. 1. Top panel: GOES light curves from the 1.0 - 8.0 A (black line) 
and 0.5 - 4.0 A (pink line) channels. Bottom panel: RHESSI X-ray light 
curves in five energy bands. 


by IWarmuth & Mantil (1201 3h also found consistent source area 
estimates among different imaging algorithms. 

Detectors 3P - 8P were employed to create images of the 
flare over the energy range 25 - 100 keV for a time period of 
approximately three minutes around the peak in HXR emission, 
23:19 to 23:22 UT. The pixel sizes were 2" x 2" and the imaging 
intervals were selected to correspond with the spectral fit inter¬ 
vals. Constructed images are shown in Pig. |2]for four intervals 
during the impulsive phase. Prom 23:20:20 UT there are two 
bright footpoints present until 23:21:00 UT, when several bright 
sources are apparent in the images. After 23:21:20 UT there is a 
sharp reduction in counts above 25 keV. The sources present in 
images after this time are broken up and diffuse. 

The HXR source area in each image was estimated as that 
contained within a certain intensity threshold. In this case, we 
adopted 30% of the highest intensity in each image. The adopted 
criteria for determining the source area provided estimates that 
ranged from 8 X lO'^ - 2x 10'^ cm^ over the time period 23:19:58 
to 23:21:20 UT. Por a single interval (23:20:44 - 23:20:56) the 
area was also estimated by fitting 2D Gaussian profiles to the 
two footpoint sources visible in the image and using the equa¬ 
tion A - nab, where a and b are the major and minor axes of 
the 2D Gaussian. The area derived using this method was found 
to be 2.4 X 10*^ cm^ for the northern footpoint and 2.2 x 10'^ 
cm^ for the southern footpoint, approximately half of the value 
estimated by summing the pixel area within the intensity thresh¬ 
old. This method was not adopted as the sources in each image 
could not always be adequately represented by a 2D Gaussian 
profile. The affect of alb edo may also have co ntributed to the 
estimated footpoint area. iBattaglia et alJ (1201 ih estimated that 
albedo could result in increased source areas by up to 30% for 
energies between 10 and 100 keV. Given these considerations, 
the measured footpoint area from HXR imaging should be con¬ 
sidered as an overestimate. 

2.2. White-iight imaging 

In addition to the source area estimates from HXR images, con¬ 
tinuum intensity obser vations made by the Helioseismic and 
Magnetic Imager (HMI, IScherrer et al.ll201^ were examined in 
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Fig. 2. Top row: RHESSI imaging of the flare on 2011 March 09 over the time range UT 23:20:20 - 23:21:08. Images were created over the energy 
range 25.0 - 100.0 keV using the Pixon algorithm with an integration time of 12 seconds and detectors 3F-8F. The applied colour scaling is based 
on the highest and lowest data values present in the images displayed. Overplotted in grey are the 10% and 30% intensity contours. Bottom row: 
Spectral fits in count flux units made during the same intervals.The background-subtracted data are plotted in black with the background level 
plotted in pink. The total fit is plotted in red along with the thermal component (green), the thick-target component (orange), albedo component 
(purple), and pulse pileup correction (cyan). The normalised residuals are shown in each plot in units of sigma. 


an attempt to identify and estimate the size of the optical flare 
kernels. A mean image was calculated from 23:15:29 - 23:18:29 
UT and subtracted from the images taken during the impulsive 
phase of the flare. Two optical kernels were identified within the 
area of the HXR footpoints. Figure [3] displays the HMl differ¬ 
ence image at 23:20:44 UT, with the 30% HXR contours over¬ 
laid in red and the 30% and 50% WL intensity contours overlaid 
in green and blue, respectively. A 10" x 10" region was centred 
around each footpoint and the area of the WL emission estimated 
as that contained within the 50% intensity contours in these re¬ 
gions. This provided estimates of the WL emitting area to be 
7.8 X 10^® to 1.0 X 10'^ cm2 between 23:20 and 23:22 UT, in 
agreement with previous studies of flare kernels in optical obser- 
vations that obtained area estimates i n the range 10'® to lO'^ cm^ 
(iKrucker et al.ll201l1: IXu et alj|2012h and an order of magnitude 
less than that obtained from HXR imaging. 

2.3. Spectroscopy 

The HXR spectroscopic analysis was performed using disk- 
integrated RHESSI spectra in Object Spectral Analysis Execu¬ 
tive (OSPEX). Spectral fits were made individually to six front 
detectors (1, 4, 5, 6, 8, and 9). The process of fitting to each 
detector individually was used to verify that the parameters re¬ 
turned from each detector were consistent with the others. The 
use of a single detector or combined data from a subset of detec¬ 
tors could have introduced a systematic error into the derived pa¬ 
rameters of the thermal and non-thermal emissions. This process 
led to detector 3 being excluded from the analysis as the results 
obtained were not consistent with the other detectors analysed. 
An average value of the parameters derived from the six detec¬ 
tors at each interval was adopted, with the error estimated as the 
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Fig. 3. HMI difference image at 23:20:44 UT. Overlaid are the 30% 
contours (red) from the closest HXR image at 23:20:50 UT and the 
30% and 50% white-light intensity contours (green and blue). 

standard deviation (iMilligan et al.ll^OMh . Spectra were fit in 12- 
second intervals from approximately 23:20 to 23:21:30 UT with 
a longer interval on either side to account for the lower count 
rate. The intervals were chosen to ensure the spectra were not fit 
across a change of attenuator state and to match with the imag¬ 
ing intervals. Spectral fits for four intervals are shown in Fig.|2] 
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Standard corrections for decimation and pu lse pileup were ap- 
plied, and the affect of photospheric albedo (iKontar et al.ll200^ 
was corrected for assuming isotropic emission. An isothermal 
and thick-target bremsstrahlung component were forward-fit to 
the spectra. The non-thermal emissions were well represented by 
a single power law over the duration of the flare. Therefore, the 
electron distribution was characterised by three free parameters: 
the spectral index {6), the low-energy cutoff {Ec) and the electron 
flux (Fq), and the thermal component by two free parameters, 
the temperature (T) and emission measure {EM). Throughout 
the duration of the flare, the spectra would be characterised as 
soft, with values of the spectral index in the range 7-8, and the 
highest emission detected above the background was typically 
70 - 80 keV around the flare peak. As a result of the very soft 
spectra, there was no apparent break between the thermal and 
non-thermal emission. 

2.4. Parameter evolution 

The evolution of the parameters derived from HXR spectroscopy 
and imaging are shown in Fig.|4l There was little variation in the 
value of the low-energy cutoff, and it remained between 22 - 
28 keV for the time period studied. The spectral index was soft 
over the entire flare duration, ranging between 7 and 8, but did 
become harder at approximately 23:20:30 UT and 23:21:15 UT 
corresponding to the two peaks in the 50 - 100 keV light curve. 
Estimates of the HXR source area reached a minimum at the 
flare peak and then increased with time as the sources became 
more extended. After 23:21:20 UT, the spectral index steepened 
and the low-energy cutoff decreased, resulting in the calculated 
energy content in non-thermal electrons increasing while at the 
same time the observed counts above 25 keV decreased. The 
parameters returned from the spectral fits made to individual de¬ 
tectors became discrepant after this time. 

Dividing the power contained in non-thermal electrons above 
the low-energy cutoff by the HXR source area provided an es¬ 
timate of the energy deposition rate into the atmosphere. The 
flux increased from 10*° erg cm“^ s“' at 23:19:30 UT (f = 0), 
reaching 10*' erg cm“^ s“' after approximately a further 50 
seconds, while the highest energy deposition rate of 2.5 x lO" 
erg cm“^ s ' occurred at 23:20:40 UT. While the WL source ar¬ 
eas indicated that the energy deposition rate into the atmosphere 
could be greater than 10*^ erg cm“^ s“*, these sources would be 
formed at photospheric heights. If the magnetic field lines con¬ 
verge towards the photosphere, then it may not be appropriate 
to use them to derive the energy deposition rate of the lower 
energy electrons stopped higher in the atmosphere. The derived 
electron parameters were adopted as input to the simulation until 
23:21:20 UT. After this time, there did not appear to be compact 
footpoint sources present in the images. 

3. Model 

The derived electron beam parameters were used as input to 
the one-dimensional NLTE radiative transfer and hydrodynam¬ 
ics modelling code RADYN. In the code, the plane-parallel 
equations of hydrodynamics are solved along with the statisti¬ 
cal equilibrium, population con servation, and radiati ve transfer 
equations on an adaptive grid (iDorfi & DrurvI [198^ that con¬ 
sists of 191 grid points with a base in the photosphere and an 
apex in the corona at a height of 10 Mm. Solutions are obtained 
for a six-level hydrogen atom, a nine-level helium atom, a six- 
level singly ionised calcium atom, and a four-level singly ionised 
magnesium atom. In total, 24 bound-bound and 22 bound-free 



Fig. 4. Evolution of the electron beam parameters determined from 
HXR spectral analysis. The quantities plotted from top to bottom are 
the emission measure, temperature, spectral index, low energy cut-off, 
total electron flux, estimated source area, power in non-thermal elec¬ 
trons above the cutoff, and the derived energy deposition rate into the 
atmosphere. The parameters from detector 1 are plotted as a solid line 
(red), detector 4 as a dotted line (black), detector 5 as a dashed line 
(green), detector 6 as a dash-dot-dash line (blue), detector 8 as a dash- 
dotted line (yellow), and detector 9 as a long dashed line (brown). 


transitions are solved in detail, including the first four transi¬ 
tions of the Lyman series, the first three Balmer transitions, the 
Ca II H & K lines and infrared triplet, the Mg ii h & k lines, and 
the He ii 304 A line. Complete frequency re-distribution was as¬ 
sumed for all transitions except the Lyman transitions, for which 
the effects of partial re-distribution are mi micked by truncat¬ 
ing th e line profiles at ten Doppler widths (iMilkev & MihalasI 
Il973h . The equations are solved using five angle points and up 
to 100 frequency points for each transition. Other atomic species 
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Fig. 5. Structure and evolution of the simulated solar atmosphere. The top panels display the atmosphere in the early phase when the beam-heating 
rate is low. The bottom panels show the evolution of the atmosphere during the explosive phase. The temperature (black) and electron density 
(red) are shown in the top row. The second row shows the pressure (black) and mass density (blue). The third row displays the beam-heating rate 
(black) and ionisation fractions of He II (green). The bulk velocity of the atmosphere (black) and ionisation fraction of He III (orange) are shown 
in the bottom row. The original atmosphere is overplotted as dotted lines at T= 100s to illustrate the change in atmospheric structure. 
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are included in the calculations as bac kground con t inua in LTE 
using the Uppsala opacity package of iGustafssonI (Il973h . Op¬ 
tically thin cooling from bremsstrahlung and coronal metals is 
included using emissivities from the CHIANTI atomic database 
dPere et al, 1_997^ |Lanxh et alj^ 201 3h . The beam-heating func¬ 
tion of Abbett & Hawley ( 19991) b^ed on the analytical treat¬ 
ment of lEmsliJil 19781) was used for an electron distribution de¬ 
scribed by a single power law. Non-thermal collisional excitation 
and ionisat i on wa s included for hydrogen using the treatment of 
lEang et akl (Il993l) . A complete description of t he method of so- 
lution can be found in fAbbett & Hawl^ (Il999h and lAllred et alJ 
(l2005h . The simulation was performed on 12 processor cores, 
and the model presented here took approximately one month 
to complete the calculations. Electron beam heating ended af¬ 
ter 110 seconds in the simulation, then the model was allowed to 
continue evolving for a further 300 seconds. 


3.1. Atmosphere structure and evolution 

The evolution of the model is plotted in Eig. |5] displaying the 
loop temperature, electron density, gas pressure, mass density, 
beam heating rate, velocity, and ionisation fractions of He ii and 
He III as a function of height in the atmosphere at ten times dur¬ 
ing the beam-heating phase of the simulation. The initial pre¬ 
flare atmosphere is structured with a chromosphere from 1 - 
1.6 Mm above the photosphere (defined as the height where 
Tsoonm = 1), the transition region (TR) at 1.6 Mm, and a corona 
with a temperature of 1 MK from 2-10 Mm (Eig. |5] top left 
panel). The location of highest non-thermal electron energy de¬ 
position was approximately 1.0 Mm above the photosphere at 
t = Os, at the base of the pre-flare chromosphere. The energy de¬ 
position rate for the first 30s was approximately lO'® ergs cm“^ 
s * and the initial phase of the flare was computed quickly as the 
atmosphere remained in the regime of gentle evaporation, with 
the majority of the energy input by the beam being balanced by 
optically thin radiative losses. This initial heating resulted in an 
upward motion of the transition region and an extended region 
of increased temperature and electron density between 1-3 Mm. 

The energy deposition rate into the atmosphere began to in¬ 
crease after 23:20:00 UT (f = 30 s) and reached 10" erg cm“^ 
s * at approximately 23:20:10 UT, at which time the explosive 
evaporation phase began. In the region of peak beam heating, 
helium had become fully ionised by f = 40s, the temperature of 
the plasma began to rapidly increase and a region of high pres¬ 
sure formed. Material expanded outwards from this region as a 
result of the pressure gradient between the cool material below 
and the low-density coronal material above (lEisher et al.lll9853) . 
A downward-directed compression wave moved through the at¬ 
mosphere with a locally cool region referred to as the “chro¬ 
mospheric condensation” (lEisher et al.lll985ah located between 
the compression wave and the flare TR. The downward-directed 
wave propagated with velocities of a few lO’s km s“’, with de¬ 
creasing velocity as it travelled into regions of higher density, 
while the explosively evaporated material propagated with ve¬ 
locities of several hundred km s '. By 45 s into the simulation, 
or approximately 5 s after the explosive phase began, the temper¬ 
ature of the lower atmosphere between 1.5-2 Mm had increased 
to 10 MK and the electron density in this region was lO" cm“^. 
This is illustrated in the bottom panels of Eig. |5]from T = 40 to 
50 s. The presence of high-temperature dense plasma located at 
low altitudes in the simulation agrees with several observational 
studies of the temperature and density of f ootpoint plasma made 
dur ing the impulsive ph ase, for instance [Hudson et akl (Il994h 
and lciraham et al.l (l2013h . 


The expansion of the hot material outwards resulted in a very 
narrow, dense region moving upward through the atmosphere. 
The very high electron densities (n^ = 10*"' cm“^) lead to a local 
temperature minimum due to the strong radiative losses from this 
material. However, after approximately 48 s of the simulation 
this region began to rapidly increase in temperature, from 10^ K 
to over 10^ K in approximately 0.1s. This is illustrated in Eig.|6] 
which displays the ion fraction of He iii, the temperature, and en¬ 
ergy balance of the atmosphere. Examining the energy balance 
of the atmosphere showed that there was a strong conductive en¬ 
ergy flux (Eig. |6] orange lines) into the cool material as a result 
of the steep temperature gradient. Radiative losses from optically 
thin emission (cyan lines) and from transitions solved in radia¬ 
tive transfer (blue lines) were insufficient to offset the energy 
input into this region, resulting in a net energy gain (black lines). 
This scenario appears similar to the initial explosive evaporation 
scenario that occurs in the chromosphere, where a high-energy 
flux leads to runaway ionisation and heating of the plasma. To 
our knowledge, this scenario has not been previously observed 
in RHD flare simulations, but the heating of this locally cool re- 
gion has been observe d in hydrodynamic flare simulations (e.g. 
iNagai & Emsliel[l984h . The very soft spectral indices and low 
values of Ec that have been employed may be responsible for this 
because a greater proportion of non-thermal energy would be de¬ 
posited higher in the corona compared to previous RHD simula¬ 
tions, which typically used harder electron distributions. As the 
temperature of the atmosphere in this cool region increased, the 
problem of solving the radiative transfer in the atmosphere was 
simplified as the transitions solved in radiative transfer were no 
longer being formed in two distinct regions of the atmosphere. 
Approximately 1.6 x 10® time steps were required to reach 48.1s 
of simulated solar time, while only a further 1.0 x 10® time steps 
were required to advance the simulation to 400s. The evolution 
of the atmosphere could then be studied over longer timescales 
compared to previous RHD flare simulations that were subjected 
to high-energy fluxes. We note that a similar atmosphere evolu¬ 
tion has been observed in subsequent simulations involving very 
soft spectral indices and/or low values of Ec. A thorough inves¬ 
tigation into the conditions that lead to this, involving a grid of 
models with varying spectral index, low-energy cutoff, and en¬ 
ergy deposition rate, is currently being undertaken. 

As the simulation progressed, the peak of non-thermal elec¬ 
tron energy deposition moved to higher altitudes as the coronal 
density increased and by 75s, the majority of energy was being 
deposited above 6 Mm. The deposition of energy in the corona 
produced a large increase in temperature. By the end of beam 
heating, the coronal temperature had increased from 1 to 50 MK 
and the mass and electron density had increased by over three 
orders of magnitude. The TR moved downward in altitude and 
the lower atmosphere was compressed. This is illustrated in the 
upper panel of Eig. |7] which displays the temperature and den¬ 
sity structure of the lower atmosphere between 0-1 Mm at 110 
s. At this time, the flare TR was located at an altitude of 500 
km above the photosphere. The temperature of the upper photo¬ 
sphere between 0.1- 0.4 Mm increased by several 100 to 1000 K 
over the course of the simulation. By examining the contribution 
functions of the transitions solved in radiative transfer, we could 
investigate the formation heights of the differe nt emission mech¬ 
anism s. The intensity contribution function dCarlsson & SteinI 
Il997h is given by Eq. [1] and represents the fractional contribu¬ 
tion to the emergent intensity originating from each height in the 
atmosphere. The terms in this equation are the source function 
S y, the monochromatic opacity per unit volume Xv , and the op¬ 
tical depth Ty specified at frequency v and a viewing angle p. 
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Fig. 6. Top panels: The ionisation fraction of He III from 48.07 to 48.14 s. Time is increasing from the left to right side of the plots. The middle 
panels display the temperature of the atmosphere. The bottom panels display the energy balance of the atmosphere. The quantities displayed are the 
net energy gain (black), viscous heating (red), work done by pressure (green), radiative losses from transitions solved in radiative transfer (blue), 
optically thin losses (cyan), conduction (orange), the background-heating function (brown), and the energy input from beam electrons (grey). 


Ci(1) 

Ty 

The bound-bound transitions were formed in an extremely nar¬ 
row (<1 km) and high-density (n^ a; 10'^ cm"^) region at the 
base of flare TR, illustrated by the sharp peak in the H a contri¬ 
bution function in the bottom panel of Fig. |7] The emitting loca¬ 
tion of these transitions changed as a function of time through the 
simulation as the TR moved downward through the atmosphere. 
The contribution function for the Balmer continuum at 3646 A is 
also plotted in this figure, illustrating that the optical continuum 
emission mainly originated from the heated layers of the upper 
photosphere with a local maximum across the flare TR. 

Once energy deposition by non-thermal electrons ceased, the 
atmosphere began to cool and material drained from the corona. 
The energy balance of the atmosphere showed that conductive 
cooling dominated the coronal energy losses for the first 50 s 
following the end of beam heating, after which radiative cooling 
dominated the coronal energy loss. After 300s of simulated time, 
or approximately 200s after the end of beam heating, the coronal 
plasma had cooled down to 1 MK and then became radiatively 
unstable. The material present in the corona rapidly decreased in 
temperature, resulting in a cool (T » 10^ K), dense (rie ~ 10" 
cm“^^) loop. 

3.2. Electron collisional stopping depth 

The changes in the mass structure of the atmosphere during 
the flare affect where non-thermal electrons will be collision- 
ally stopped and lose their energy. We investigated where in the 
atmosphere electrons would be stopped at different times in the 



Height (Mm) 

Fig. 7. Top panel: Structure of the lower atmopshere at 110s into the 
flare simulation. Solid black and red lines display the temperature and 
electron density. The dashed lines display the original atmosphere at t 
= Os. Bottom panel: contribution function for the Balmer continuum at 
3646 A and for H a at line centre (6564.7 A) for a viewing angle of 
p = 1. The H a function was scaled by a factor of 100 to be visible 
on the same scale. The dashed black line displays the inital contributon 
function for the continuum at t = 0 s. 
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evolution of the simulated flare. The collisional stopping column 
density of electrons with energies fro m 5 to 200 keV was calcu¬ 
lated using the analytical treatment of lEmslij(Il978h (Eq.|2]l, and 
the corresponding height in the atmosphere was then determined 
from the flare simulation. The results are shown in Eig. 0 for 
three times in the simulation at f = 0, 70 and 110 s, 

I™'’- ® 

The black line in Eig. 0 shows that the stopping height in the 
initial atmosphere for the energy range considered lies between 
0.8 to 1.5 Mm above the photosphere. As the density of the 
corona increases, lower energy electrons are stopped higher in 
the atmosphere, with the majority of electrons being stopped 
near the apex of the coronal loop. If we examine the stopping 
depth of electrons of a specific energy, for example 50 keV, we 
can show how the greatest penetration depth changes through¬ 
out the simulation. In the initial atmosphere, 50 keV electrons 
would be stopped at approximately 1.1 Mm above the photo¬ 
sphere, while at 70s they can penetrate to around 0.75 Mm, as 
the column density required to stop the 50 keV electrons is now 
located at a lower altitude in the atmosphere. At the end of beam 
heating, the 50 keV electrons would now be able to penetrate 
no farther than 5 Mm into the atmosphere, as the density of the 
corona has increased. This illustrates that the column density re¬ 
quired to collisionally stop an electron corresponds to different 
altitudes above the photosphere at different times in the simula¬ 
tion. The greatest possible penetration depth was approximately 
0.45 Mm above the photosphere at 130s into the simulation for 
electrons of energy 70 keV and higher. They were stopped at the 
base of the flare TR where there is a sharp increase in the col¬ 
umn density to over 10^^ cm“^. These low altitudes are similar 
to WE observations made in the contin uum at 6173 A of a limb 
flare by iMartmez Oliveros et al.l (l2012h . who found WL source 
heights of 0.3 Mm above the photosphere and co-spatial with 
30 - 80 keV HXR emission. These source heights were located 
below the expected penetration depth of the 50 keV electrons 
that would be required to produce the observed HXR emission. 
The downward mass motion and compression of the lower at¬ 
mosphere in our simulation could offer an explanation for the 
co-spatial WL and HXR emission present at low altitudes in this 
observed flare. 

Erom the analytical formula for the electron collisional stop¬ 
ping column density we have demonstrated that it is possible 
for energetic electrons to propagate farther into a flaring atmo¬ 
sphere compared to a quiet-Sun atmosphere. This is explained 
by material being compressed below the altitude at which explo¬ 
sive evaporation occurs. In general, it was found that electrons 
with energies below 50 keV are collisionally stopped at a higher 
altitude as material is ablated into the corona, while the more en¬ 
ergetic electrons (>60 keV) are able to propagate farther into the 
atmosphere. It required continuous energy input into the atmo¬ 
sphere over a timescale of 100 s to produce these deep penetra¬ 
tion depths. The simulation presented here used very soft spec¬ 
tral indices and moderate energy deposition rates. The use of 
higher energy deposition rates may produce a more rapid com¬ 
pression of the atmosphere, allowing electrons to penetrate to 
low altitudes on shorter timescales. Einally, we note that this 
analysis is quite simplified, using the analytical expressions for 
electrons penetrating into a fully ionised cold-target. A more ad¬ 
vanced treatment using the Eokker-Planck version of the RA- 
DYN code may alter the depth at which the electrons are able to 
penetrate to, as it includes pitch-angle scattering and magnetic 
mirroring. 



Electron Energy (keV) 


Fig. 8. Greatest penetration depth of electrons injected into the atmo¬ 
sphere as a function of energy from 5 - 200 keV for three times in the 
simulation. The inset figure displays a zoom-in of the atmosphere be¬ 
tween 0-2 Mm. 


4. Observational comparison 

The results of the model, including hydrodynamic variables and 
radiative emissions, were compared to observations made with 
EVE, GOES, and RHESSI. A pre-flare background averaged 
over a five-minute period was subtracted from the EVE data to 
obtain the flare irradiance. 


4. 1. Hydrodynamic variables 


The evolution of the hydrodynamic variables during the studied 
flare were measured using spectroscopic observations from the 
EVE MEGS-A instrument. The Ee xxi 145/128 density-sensitive 
line ratio w as used to determine the mean electron density of 
the plasma dMilligan et al.ll^l2bh . Emission lines from Eeviii 
to Eexxiv allowed for the temperature structure of the flaring 
plasma to be est imated from a differ ential emission measure 
(DEM) analysis (iKennedv et al.ll2013). The DEMs were con¬ 
structed using coronal abundances dEeldman et al.lll992h over 
the temperature range log T^ = 6.0 - 7.5. Estimates of the 
isothermal T and EM were obtained from GOES filter ratios 
(IWhite et alJl2005h available through the goes widget in Solar- 
Soft, and the T and EM were also determined from analysis of 
RHESSI spectra. The evolution of the electron density, T and 
EM were measured and compared to the modelled parameters. 
Detailed descriptions of the diagnostics methods can be found in 
the papers referenced above. The model coronal electron density 
was calculated in a manner to directly compare with the mean 
electron density determined from Eexxi density-sensitive line 
ratios. An average value was calculated using the grid cells that 
fell within the typical formation temperature range of Ee xxi (log 
Tg = 6.8 - 7.2). The value of tig at each grid cell was weighted 
by the height Az of that cell, to avoid the average density being 
dominated by regions of the highest spatial resolution. 

Eigure 0 displays the modelled and observed values of the 
EM, T, and electron density. The temperature measurement from 
EVE is the peak temperature of the DEM at each time bin, and 
the EM is the integral of the DEM over the range log Tg = 6.0 - 
7.5. Values of the flare emission measure obtained from RHESSI 
and GOES are very similar. However, the EM derived from EVE 
is consistently lower than both of the RHESSI and GOES mea¬ 
surements. The model EM was the integral of the model DEM 
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above log > 6.0. It was scaled by an estimate of the emitting 
area and has a peak value of log EM = 49.1 after 150 s into the 
simulation. The peak temperature determined from all three in¬ 
struments was approximately log = 7.3 (20 MK), while the 
peak loop temperature of the model reached log Tg = 7.7 (50 
MK). The combined EVE and RHESSI DEM for this event was 
studied bv ICasni et alJ (l2014l) . who also found peak temperatures 
between 20 - 30 MK with no evidence for a super-hot compo¬ 
nent. These observations indicate that the coronal temperature in 
the model is too high. A possible explanation for this is that the 
value of Ec determined from the spectral fitting was too low, re¬ 
sulting in a large amount of energy being contained in lower en¬ 
ergy electrons that deposit their energy higher in the atmosphere. 
Alternatively, it might mean that the continuous deposition of 
energy into a single model atmosphere for a prolonged period of 
time is unrealistic. The comparison with the observed electron 
density was limited as the density-sensitive lines were weak in 
the EVE spectra resulting in large irradiance uncertainties, and 
were also affected by line blending. The modelled Fe xxi electron 
density reached a peak value of log — 11.8 after 100s of the 
simulation followed by a slow decay over several minutes. EVE 
measurements indicated values of log n^ = 11.0 - 11.5, which 
were slightly lower than the model atmosphere. However, the 
density at the apex of the modelled loop (see Fig. 0 agreed bet¬ 
ter with the observed electron density, reaching a peak value of 
approximately log n^ = 11.5. 

The temporal evolution of the modelled parameters was also 
quite different from what was observed. Pre-impulsive phase 
heating leads to high-temperature plasma being present in the 
corona before a detectable increase in HXR emission, and the 
temperature and EM decrease much more rapidly than in the 
observations. Additional energy release in the decay phase that 
has not been accounted for in our simulation could explain this. 
This would agree with previous numerical modelling and obser¬ 
vational studies of decay phase cooling, which indicated that ad¬ 
ditional energy release is required to explain the cooli ng times of 
flare lo ops, for instance fReeves & WarrenI (l2002h and iRvan et alJ 
(l2013h . It is possible to add a time-dependent volumetric heating 
rate to the simulation to account for additional energy release in 
the decay phase of the flare. This would lead to a more gradual 
decay as the additional heating would offset the radiative losses 
from the atmosphere. However, we have not attempted to model 
this without a physical description of the additional energy de¬ 
position in the atmosphere. 

4.2. Light curves and synthetic spectra 

The only bound-bound transition that is solved in NLTE radia¬ 
tive transfer in RADYN and was observed during this event is the 
He II 304 A doublet. This transition and the He ii recombination 
continuum are observed by EVE MEGS-A. The flux of the 304 
A doublet was measured by fitting with a Gaussian profile, while 
the flux of the continuu m was measured using a method simi¬ 
lar to that described by iMilligan et al.l (1201 2ah . by fitting with 
an exponential function and subtracting the underlying free-free 
emission. Light curves of EUV emission lines and the He ii dou¬ 
blet were compared to synthetic light curves generated from the 
RADYN model. The RADYN light curves are displayed using 
all output time steps, but to facilitate direct comparisons with 
observations, the light curves are also displayed degraded to the 
cadence of the observing instrument. Flux was averaged over the 
instrument exposure time and sampled at the instrument cadence 
(both 10 seconds for EVE). 



23:20 23:21 23:22 23:23 23:24 

StartTime(09-Mar-ll 23:19:30) 


Fig. 9. Comparison of the modelled and observed emission measure, 
temperature, and electron density. The red lines are the measurements 
from GOES, blue lines from EVE, brown lines from RHESSI, and black 
lines are the modelled values. 

The modelled and observed light curves of the He ii recom¬ 
bination continuum and He ii 304 A are shown in Fig. [TO] The 
degraded light curve is plotted in blue. The modelled values are 
for a viewing angle of p - 0.95 and have been scaled by an 
estimate of the flare area then converted into flux that would 
be measured at Earth to facilitate comparison with the observa¬ 
tions. The observed light curves are displayed below each of the 
modelled light curves. In the hrst 45s of the simulation the light 
curves feature sharp, sub-second variations in intensity. These 
variations are smoothed when the light curve is degraded to the 
instrumental exposure time. After the cool material moving up¬ 
ward through the corona increased in temperature, the intensity 
of He II declined to a minimum as helium became fully ionised 
in this region. The intensity then increased to a maximum as 
the density in the region of line formation at the base of the 
TR peaked at approximately 75 s. The time of peak intensity of 
the observed and modelled He ii 304 A light curves occurred to 
within 40s of each other. The modelled light curve decayed faster 
as energy input into the atmosphere ceased after 100s of the sim¬ 
ulation. The modelled 304 A emission was an order of magni¬ 
tude greater than the measured flux of the line by EVE, and the 
modelled continuum was approximately 2 orders of magnitude 
greater than the observed continuum flux. 

Synthetic EUV spectra were generated by constructing 
DEMs of the simulated atmosph ere averaged over 10 s intervals. 
The CHIANTI atomic database dLandi et al.l2013h was then em¬ 
ployed to generate spectra of the 60 - 370 A wavelength range at 
the wavelength binning (0.2 A) of the EVE MEGS-A instrument 
and at the full width half maximum of observe d lines in EVE 
spectr a (0.7 A), assuming coronal abundances (iFeldman et al.l 
Il992h . The comparison between the observed and synthetic EUV 
emission is displayed in Fig. [TT] at three times during the flare 
evolution. The best agreement between the modelled and ob¬ 
served EUV emission is found at approximately 100 - 140 s into 
the simulation (Fig.[TT] middle panel). At these times, the mod¬ 
elled irradiance agreed with the observed irradiance to within a 
factor of two. Considering that the calculation of the synthetic 
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Fig. 10. Light curves of He ii recombination and 304 A line emission 
created from the RADYN model compared to observed light curves 
from EVE. The modelled light curves are displayed at each saved time 
step and then degraded to the instrumental cadence and exposure time 
of EVE (plotted in blue). The RADYN light curves have been scaled by 
an estimate of the flare area and converted to the flux measured at Earth 
to facilitate direct comparison with the observations. The observed light 
curves from EVE MEGS-A are displayed at 10s cadence. 
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Fig. 11. Observed (black) and synthetic (red) spectra of the 90 to 200 
A wavelength range. The prominent emission lines are the flare lines of 
Fe XVIII to Fe xxiii from 90 - 150 A, emission lines of Fe viii to Fe x near 
165 - 180 A, and the Fexxiv line at 192 A. 


irradiance relied on the assumption of ionisation equilibrium, el¬ 
emental abundances, and an estimate of flare area, this is reason¬ 
ably good agreement. The irradiance of emission lines formed 
at temperatures of between 0.1 and 1 MK was overestimated by 
around a factor of ten during the impulsive phase and peak of 
the flare. The overestimation of the intensity of UV em ission 
lines was also found in the models of iFisher et all (Il985bl) . This 
is due to the high electron densities across the model flare TR, 
which resulted in an extremely high EM at these temperatures, 
despite the narro w grid cell heights. During the gradual ph ase 
coronal dimming (iRust & HildneijlT976t IWTOds et al.ll201lh af¬ 
fected the spatially integrated EVE observations and reduced the 
observed irradiance of these lines, preventing an accurate com¬ 
parison to the modelled irradiance. Light curves of several emis¬ 
sion lines were generated by fitting the synthetic spectra with 
Gaussian profiles in the same manner as observed EVE spec¬ 
tra are analysed. The line irradiance and emission profiles were 
compared to what was observed during the flare (see Eig. [I2]i. 
Emission from Ee xxiv line peaked first at 130s into the simula¬ 
tion with emission from cooler lines peaking later as the corona 
cooled. The Ee xxiv emission peaked approximately 100s before 


the observed light curve, with a longer delay between the time 
of observed and modelled peak flux for cooler emission lines. At 
330s into the simulation, the loop became radiatively unstable 
and cooled to approximately lO'^ K. 

5. Conclusions 

We have simulated the response of the solar atmosphere to non- 
thermal electron beam heating using a ID RHD modelling code. 
Our aim was to compare the modelled and observed flare param¬ 
eters and line irradiance and to investigate the structure of the 
lower atmosphere during the flare. During the explosive evapo¬ 
ration phase the region of cool, dense material moving upward 
through the atmosphere rapidly increased in temperature. De¬ 
spite the very high electron densities, radiative losses were un¬ 
able to balance the energy input into this region via thermal con¬ 
duction. The very soft spectral indices and low values of Ec that 
were employed could explain this, due to the greater proportion 
of non-thermal electron energy that would be deposited higher 
in the atmosphere. This greatly simplified the problem of solving 
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Stai-t Time (09-Mar-l 1 23:19:30) 


Fig. 12. Observed (black) and synthetic (red) light curves of emission 
lines from Fe ions observed by EVE. The dotted lines indicate the times 
corresponding to the spectra displayed in Eig. m 

the radiative transfer in the atmosphere and allowed the evolution 
of the atmosphere to be studied over longer timescales compared 
to previous RHD flare simulations involving very high energy 
deposition rates. 

It was found that the structure of the lower atmosphere be¬ 
came extremely compressed during the simulation. The flare TR 
moved downwards in altitude during the explosive phase and 
was located at 440 km above the photosphere at 120 s. By cal¬ 
culating the collisional stopping depth of electrons, it was deter¬ 
mined that the downward mass motion and compression of the 
lower atmosphere allowed for higher energy electrons to pene¬ 
trate to a lower altitude above the photosphere when compared 
to the pre-flare atmosphere. The location of deepest penetration 
coincides with the region where line emission originates and 
with a peak in the intensity contribution function of optical con- 
tinua. This could help to explain observations of low-altitude, co - 
spatial WL and HXR emission dMartmez Oliveros et al.ll2()T^ . 
which was considered contradictory to the standard thick-target 
model. However, continuous heating of the atmosphere over a 
timescale of 100s was required to produce these electron pene¬ 
tration depths. 

High-temperature plasma was produced at the footpoint of 
the modelled loop during the impulsive phase. Less than 10s af¬ 
ter explosive evaporation occurred the temperature in the lower 
atmosphere increased to 10 MK at altitudes of 1 - 2 Mm. The 
electron density in this region was approximately 10^' cm“^. 
This generally agrees with observational studies of the tempera¬ 
ture and density of flare ribbons and footpoints, which indicate 
the presence of high-temperature material present at low alti¬ 
tudes in the solar atmosphere. The model produced a denser and 


hotter flare loop than what was detected and the temporal evo¬ 
lution of the loop temperature, electron density, and EM did not 
agree with observations. The modelled loop cooled more rapidly 
than what was observed, suggesting additional energy release 
in the decay phase. The comparison between observed and syn¬ 
thetic EUV irradiance indicated good agreement during the peak 
of the simulated flare, to within a factor of 2 for the majority 
of emission lines. However, the irradiance of lower temperature 
emission lines from ions formed at temperatures lower than 10® 
K was significantly overestimated during the very early and the 
decay phase of the simulated flare. 

Although the energy deposition rate into the atmosphere is 
an underestimate, modelling of this flare has highlighted several 
interesting results. A grid of flare models is currently being gen¬ 
erated to investigate the atmospheric response to heating from 
different distributions of non-thermal electron energy across the 
observed range of the parameters that describe the electron dis¬ 
tribution. The inclusion of higher energy fluxes and the adoption 
of different pre-flare atmospheres should also be investigated to 
determine how this affects the chromospheric response to non- 
thermal energy input and the evolution of the flare atmosphere. 
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